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ANISOTROPIC MESH QUALITY MEASURES AND ADAPTATION 
FOR POLYGONAL MESHES* 

WEIZHANG HUANGt AND YANQIU WANG* 

Abstract. Anisotropic mesh quality measures and anisotropic mesh adaptation are studied 
for polygonal meshes. Three sets of alignment and equidistribution measures are developed, one 
based on least squares fitting, one based on generalized barycentric mapping, and the other based 
on singular value decomposition of edge matrices. Numerical tests show that all three sets of mesh 
quality measures provide good measurements for the quality of polygonal meshes under given metrics. 
Based on one of these sets of quality measures and using a moving mesh partial differential equation, 
an anisotropic adaptive polygonal mesh method is constructed for the numerical solution of second 
order elliptic equations. Numerical examples are presented to demonstrate the effectiveness of the 
method. 

Key words, anisotropic, adaptive mesh, polygonal mesh, generalized barycentric coordinates. 

AMS subject classifications. 65M60, 65M60. 


1. Introduction. During the last decade polygonal/polyhedral meshes have 
gained considerable attention in the scientific computing community partially due 
to their flexibility to deal with complicated geometry, curved boundaries, and local 
mesh refining and coarsening and also to their connection with the Voronoi diagram. 
Various methods have been developed for the numerical solution of partial differential 
equations (PDEs) on polygonal/polyhedral meshes, including the mimetic finite dif¬ 
ference method (see the recent survey paper ro. conforming finite element methods 
using generalized barycentric coordinates mmmmmmm, the finite volume 
method [19] [42], the virtual element method (see [45] and references therein), the dis¬ 
continuous Galerkin method Sim Eg, and the weak Galerkin method [MlHZ]. For 
the error analysis of these methods, shape regularity of polygonal/polyhedral meshes 
has also been studied and a number of shape regularity conditions have been proposed 
The minimal requirement of these shape regularity conditions unan¬ 
imously states that, for the isotropic case, each polygon or polyhedron cannot be too 
“thin”' 

In practice, one typically wants the mesh elements to be as regular in shape and 
uniform in size as possible under certain metrics, since the approximation error and 
the condition number of the discrete system depend closely on the mesh geometry. For 
simplicial meshes, this dependence has been well studied under the topic of “mesh 
quality measures” and many different computable quality measures have been de¬ 
signed; e.g., see 0 eh ed Eg. More importantly, mesh quality measures often 
play an essential role in constructing efficient mesh adaptation algorithms since they 
also provide a definition of optimal meshes. 

Traditionally, isotropic mesh adaptation has been studied, where the shape of 
mesh elements is kept as regular (in the Euclidean metric) as possible. On the other 


* W.H. was supported in part by the NSF under Grant DMS-1115118 and the University of Kansas 
GRF Award # 2301056. Y.W. was supported by the FY2015 Spring Travel Program from the College 
of Art &; Science, Oklahoma State University, during her visit to the University of Kansas in June 
2015. 

tDepartment of Mathematics, The University of Kansas, Lawrence, KS 66045. E-mail: 
whuang@ku.edu 

* Department of Mathematics, Oklahoma State University, Stillwater, OK 74078. E-mail: yan- 
qiu.wang@okstate.edu 


1 



hand, anisotropic mesh adaptation allows elements to have large aspect ratio. The 
key is to keep elements to be aligned in some extent with the geometry of the solution. 
An anisotropic mesh can be viewed as a uniform one in a Riemannian metric, as in 
the M-uniform mesh approach of anisotropic mesh generation [23j . It has been amply 
demonstrated that a properly generated anisotropic mesh can be much better aligned 
with the geometry of the solution and lead to a much more accurate solution than 
an isotropic mesh; e.g., see among many other work. 

It is worth pointing out that all of these work is concerned with simplicial meshes. 
There do not seem to exist any systematic studies on anisotropic quality measures 
and adaptation for polygonal/polyhedral meshes. 

The objective of this work is to develop anisotropic mesh quality measures for 
polygonal meshes and an anisotropic adaptive polygonal mesh method. An efficient 
and commonly used way to generate a polygonal (or polyhedral) mesh is through 
the Voronoi diagram. The latter is a partition of a spatial domain consisting of 
polygonal/polyhedral cells associated with a given set of points called “generators” or 
“seeds” such that each cell contains only points closer to its associated generator than 
to the others. Generally speaking, the quality of Voronoi meshes is not necessarily 
good. The most common strategy to obtain a high quality polygonal mesh is to start 
with a Voronoi mesh and then employ algorithms such as Lloyd’s algorithm to obtain a 
centroidal Voronoi tessellation (CVT) jTDj. CVT can also be generated under a scalar 
metric function defining distance among points, i.e., a scalar weight function specifying 
the mesh density across the domain. This results in isotropic polygonal meshes, while 
our work here focuses on measuring and generating anisotropic polygonal meshes of 
high quality. 

A major difference between triangular and polygonal meshes is that, all of the 
elements in a triangular mesh are affine similar to a single reference triangle and their 
quality in shape and size can be measured by comparing them to the reference ele¬ 
ment. Unfortunately, this does not work for polygonal meshes since elements with the 
same number of vertices cannot be mapped to a single element of the same number 
of vertices under affine mappings. The difference poses great difficulty in studying 
polygonal mesh quality measures, as one can no longer use many techniques devel¬ 
oped for triangular meshes. In this work, in parallel to the so-called alignment (for 
regularity) and equidistribution (for size) measures for simplicial meshes (50 [2], we 
shall develop three sets of polygonal mesh alignment and equidistribution measures: 
one based on least squares fitting, one based on generalized barycentric mapping, 
and the other based on singular value decomposition of edge matrices. Numerical 
tests show that all three sets of mesh quality measures give good indicators of the 
actual anisotropic mesh quality under given metrics, with individual emphases on 
slightly different aspects. Based on one of these sets of quality measures and using 
the so-called Moving Mesh PDE (MMPDE) moving mesh method [26], we construct 
an anisotropic adaptive polygonal mesh method for the numerical solution of second 
order elliptic equations, with the aim of minimizing the Lr norm of approximation 
errors. Numerical examples will be presented to demonstrate the effectiveness of the 
method. 

The paper is organized as follows. Section [2] is devoted to the development of 
three sets of alignment and equidistribution quality measures for polygonal meshes. 
An anisotropic adaptive polygonal mesh method is constructed in Section [3] based on 
one set of the quality measures and the MMPDE moving mesh method. Numerical 
results obtained with the proposed method for two examples are presented in Section 
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21 Conclusions are drawn in the last section. 

Throughout the paper, the short name “n-gon” stands for a polygon with n 
vertices. We consider meshes consisting of convex polygons in this work. We also 
assume that all polygons in a mesh are non-degenerate, i.e., any three vertices of a 
polygon do not lie on the same line. 

2. Anisotropic mesh quality measures. In this section we study three sets of 
mesh quality measures for a general polygonal mesh T on a two-dimensional bounded, 
polygonal domain Q under the metric specified by a given tensor M. The metric 
tensor M = M(x) is assumed to be symmetric and uniformly positive definite on f 1. 
It depends on the physical solution in the case of mesh adaptation for the numerical 
solution of PDEs (See Section [3]). The simplest case is M = I (the Euclidean metric) 
where the mesh quality measures evaluate the shape regularity and size uniformity of 
the mesh. 

2.1. Continuous-level mesh quality measures. Different from triangles, n- 
gons with n > 3 in a polygonal mesh generally are not affine similar to a single 
reference n-gon. The idea of measuring the quality of a mesh through comparing its 
elements to a single reference element does not work for polygonal meshes. Neverthe¬ 
less, the idea can be used indirectly. In the following, we shall take a close look at 
the development of mesh quality measures for a triangular mesh and then establish 
continuous-level mesh quality measures which will serve as a unified framework for 
developing mesh quality measures for a polygonal mesh. 

We assume that a computational mesh Tc has been chosen. We emphasize that, 
in addition to a conventional mesh, Tc can be chosen as a collection of polygons. 
The only requirement is that, for any T £ T, there exists a corresponding reference 
polygon Tc £ Tc with an associated bijection Ft : Tc —> T. The global bijection F 
is the collection of all local Ft’s , and we want it to be piecewise differentiable so that 
its Jacobian matrix exists almost everywhere. Generally speaking, Tc should consist 
of good quality polygons under the Euclidean metric. 

We now consider the mesh quality measures for a triangular mesh. To this end, 
we assume that both T and Tc are triangular meshes. Then, one can simply set Ft ■ 
Tc —> T to be an affine mapping and evaluate the quality of T with reference to Tc- 
Specifically, the triangles in T are said to have a good shape if they, when measured 
in the metric M, are similar to their counterparts in Tc measured in the Euclidean 
metric. They are said to have the same size when the ratio of their area (measured 
in the metric M) to the area of their counterparts in Tc stays constant. Denote the 
vertices of T and Tc by X; and £ t , i = 1,2,3, respectively. If we approximate the 
metric tensor M by its average on T, 



( 2 . 1 ) 


where |T| is the area of T, then the similarity requirement between T and Tc can be 
expressed mathematically as 



where er-r is a constant on T. Using the affine mapping Ft '■ Tc —> T and its Jacobian 
matrix J-r (which is constant on T), we can rewrite the above condition as 



3 


( 2 . 2 ) 


It can be shown (cf. 


Lemma 4.1.1]) that the above condition implies 
I^-Mt-It- = gtI- 


This condition should apply to all triangles in T for them to have a good shape with 
reference to their counterparts in 7c- Taking determinant and square root of both 
sides of (2.2) and using the fact that det(lT) = ITl/lTc], one gets 

|T|det(M r )i 


\Tc\ 


= (if- 


Then the uniform size requirement simply implies that <jt should stay the same for 
all the triangles. By changing its subscript T to h , i.e., <tt —> dh, we can rewrite the 
above equation into 

|T|det(Mr ) 2 = dh\Tc\- 

Summarizing it over all triangles, we can find dh as 

d h = j^-T l^lv/detCMr), 


(2.3) 


Ter 


where |7c| denotes the total area of 7c- Moreover, (2.2) becomes 

= dhl, VT £ T- 


(2.4) 


From this the so-called alignment and equidistribution quality measures can be de¬ 
veloped for a triangular mesh (cf. (2.10) below). 

We now find a continuous analogue of \2A\ . To this end, we define the computa¬ 
tional domain as 


n c = |J t c . 

Tc&Tc 

Note that Cl c can be a conventional two-dimensional domain or a collection of poly¬ 


gons. It is not difficult to see that a continuous analogue of (2.4) is 


1*MI = dl, Vx £ Cl 


(2.5) 


fL 


ci. 


where d is a constant, I = ^ is the Jacobian matrix of the mapping T 
and £ and x are the coordinates of Cl c and Cl, respectively. At the continuous level, 
the mesh is represented by the mapping J- from Cl c to Cl. Moreover, the constant d 


in (2.5) is determined by compatibility. Indeed, taking determinant and square root 
on both sides of (2.5) and integrating the resultant equation over f2 c , we have 




( 2 . 6 ) 


where fi c denotes the area of Cl c 


The condition (2.51 implies that all of the eigenvalues of matrix J 4 MJ are equal 


to d. The requirement is equivalent to the combination of two conditions: all of the 
eigenvalues are equal to each other and the square root of the product of the eigenval¬ 
ues is equal to d. These two conditions give rise to the alignment and equidistribution 
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conditions [21] . From the inequality of the arithmetic and geometric means it is not 
difficult to show that the alignment condition at the continuous level, 

^trace(j‘MJ) = dett^MJ) 1 / 2 , (2.7) 

is equivalent to requiring all of the eigenvalues of J‘MJ to be equal to each other. A 
discrete version of this condition for a triangle T requires that T (measured in the 
metric M-r) be similar to Tq- 

On the other hand, notice that the product of the eigenvalues of J*MJ is equal to 
the determinant of the matrix. We can express the second condition (the equidistri- 
bution condition at the continuous level ) mathematically as 

det(l)i/det(M) = a. (2.8) 


For a triangular mesh, the left-hand side is proportional to the area of a triangle under 
the metric M and the condition requires all of the triangles to have the same size. 


Using (2.7) and (2.8), we can define the measuring functions as 


9aii(x) = 


trace(1* 


2det(J t MJ) 1 / 2 ’ 


9e 9 ( X ) = 


det(J)\/det(M) 


(2.9) 


In the ideal case when perfect alignment and equidistribution are reached everywhere, 
one has q a n = 1 and q eq = 1 over the entire 0. For the general situation, we may 
not have perfect alignment or equidistribution. In this case, functions q a u{ x) and 
g eg (x) indicate how closely the alignment and equidistribution conditions are satisfied 
at point x. The overall alignment and distribution measures can then be defined as 


Qau = masg nK (x), 
x£f2 


Qeq = max q eq (x). 


( 2 . 10 ) 


From the inequality of the arithmetic and geometric means and the fact that 



Qeqd$, — |fi c |, 


it is not difficult to show that both Q a u and Q eq have the range [l,oo). Moreover, 
the smaller they are, the better the alignment and equidistribution of the mapping 
T. 


The above defined continuous level mesh quality measures will serve as a unified 
framework for defining actual mesh quality measures for a polygonal mesh. Clearly, 
Qaii and Q eq depend on the metric tensor M and the computational mesh 7c (which 
determines fl c ). Recall that 7c can be a conventional mesh or a collection of polygons 
and the corresponding Q c is a conventional domain or a collection of polygons. From 
(2.9), one can see that the determining factor for Q a u and Q eq is the Jacobian matrix 
1 of the bijection T and how f l c looks has little weight on the definition. As long as 
there is a well-defined 1 almost everywhere in f2, Q c does not have to be a continuous 
manifold. 

It is useful to point out that, the mapping T is defined through the local mapping 
Tt '■ 7c —> T for all T € T. Unlike triangular meshes, there does not exist an affine 
mapping between T and Tq in general for polygonal meshes. Thus, there are many 
ways to specify J- for a given polygonal mesh 77 Commonly, J~t is chosen to be a 
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non-affine mapping for a general Tq or an affine mapping for a specially chosen Tc; 
see the next three subsections. 

We now study three discretizations/discrete approximations to the continuous 
mesh quality measures. The key components are the choices of 7c and the corre¬ 
sponding mapping T (to represent the mesh T). 

2.2. Approximation 1: use least squares fitting and edge matrices. In 

this approximation, we assume that the computational mesh 7 ~c has been chosen to be 
a conventional mesh or a collection of reference polygons. Typically 7c should consist 
of good quality polygons in the Euclidean metric. These good quality polygons can 
be regular n-gons or almost regular n-gons. 

As part of the definition, we choose to approximate the metric tensor M on T 
by its average This approximation is reasonably accurate while making M be 

constant on each element. To have a constant M on each element is important since 
it can significantly simplify the analysis and derivation. 

As mentioned before, T is defined through the local mapping Tt '■ Tc —> T. The 
existence of such a local mapping is guaranteed by the Riemann mapping theorem, 
which states that there exists a conformal mapping between any two n-gons. As we 
will see in the next subsection, there actually exist many such mappings. They all 
satisfy 


= * = 1) 


( 2 . 11 ) 


where and x^, i = 1, • • • n, are the vertices of Tc and T, respectively. Instead 
of constructing a specific example of these mappings, we consider to fit an affine 
mapping At£ + c in the least squares sense based on the above conditions. Note that 
such a mapping is an approximation to all bijections between Tc and fl satisfying 
(2.11). Although the approximation is very rough, it serves our purpose to define 
mesh quality measures. The condition for the least squares fitting is 

Ar^i + c = Xi, i = 1, ■■■, n. 


Subtracting the first equation from the others, we get 

AriZi ~ £r) = Xj — xi, i = 2, ...,n. 

A least squares solution for the above condition is 

At = E t E^ c , 

where Et and Et c are the edge matrices defined by 

E t = [x 2 - x 1} • • • , x„ - xx] <E K 2x( ” _:L) , 

Era = [£2 -€r, •••, €n - £l] e ]R 2x( " -1) , 

and E^ c is the pseudo-inverse of Et c ■ Since both T and Tc are convex and non¬ 
degenerate, both Et and Et c have full row rank. Then, the pseudo-inverse of Et c 
can be found as 


E+ c = E^E^E ^)- 1 g K (n " 1)x2 , 

which gives rise to 

At = EtE^^EtcEt ^- 1 . 
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( 2 . 12 ) 



By construction, At is an approximation to the Jacobian matrix of any piecewise 
continuously differentiable bijection between Tc and T satisfying (2.11). 

Replacing 1 and M by At and M t in (2.9) and (2.10), we obtain the first set of 
mesh quality measures for a polygonal mesh, 


(2.13) 

(2.14) 

where 

<Th,i = jr'52 det (EtE^ ( Et c £t c r 1 ) v / det(M T ) (2.15) 

P tgT 


^ aH ’ 1 “Jr 2 det ( [E t E^ c (E Tc E t Tc ) _1 ] 4 M-r [E T E^ c (E Tc E ^ c ) _1 ]) 1 / 2 ’ 

det {EtE^ {Et c Ej, )“ 1 ) \/det(MT) 

Qeq.i = max- k - k -, 

Ter a hl 


and N p is the number of the polygons in T- Notice that Uh,i is not a direct approx¬ 
imation of a defined in (2.6). Instead, it is defined by summarizing the following 
discrete equiditribution condition over all polygons, 


det (E t Et c (E Tc Et c ) 1 )\/ det (Mp) — 


The same strategy will be used in defining other two sets of mesh quality measures 
without explicit mentioning. This also makes the current formulas consistent with 
those developed for simplicial meshes in mm- 

To show how well Q a u, l and Q e q, l work, we present numerical results obtained 
with Lloyd’s algorithm for centroidal Voronoi tessellations (CVTs) [32]. The algorithm 
is known to produce a sequence of Voronoi meshes that converges to a CVT. We expect 
to see that both Q a u,i and Q eq ,i (with M = I) for the sequence of meshes decrease 
and converge to one if they are correct indicators for the shape regularity and size 
uniformity of the Voronoi meshes. 

In Fig. |2.1[ we show Voronoi meshes of 8 x 8, 16 x 16 and 32 x 32 cells obtained 
with Lloyd’s algorithm. Visually we can see that the meshes are becoming better 
in the sense that the cells are getting more regular in shape and more uniform in 
size. To compute Q a ii, i and Q eq , i, we choose the unitary regular n -gon with vertices 
£.j = (cos^P, sinpP) 4 , i = 1, ...,n as Tc for each n-gon in T. The results are 


reported in Fig. 2.2 One can see that both Q a u, i and Q eq ,i decrease towards one 


as the number of iterations increases. This reflects the convergence nature of Lloyd’s 
algorithm. Moreover, the decrease is more significant in the first few iterations and 
then getting much slower. This confirms the fact that Lloyd’s algorithm is efficient 
for obtaining a CVT of reasonably good quality but very slow for a CVT with high 
accuracy. Furthermore, one may notice that the decrease of Q a u, i and Q eq .i is not 
monotone. This may be attributed to the facts that (a) Lloyd’s algorithm is not 
designed specifically to minimize these quality measures and (b) Q a ii,i and Q eq ,i 
measure the quality of worst mesh elements. Overall, we see that Q a u, i and Q eq , i 
correctly reflect the polygonal mesh quality under the identity metric. 


2.3. Approximation 2: use generalized barycentric mappings. This ap¬ 
proximation is similar to Approximation 1 except that here we construct a specific 
local mapping Ft ■ Tc —> T using generalized barycentric mappings m uni nn m\- 
The generalized barycentric mappings are related to generalized barycentric coordi¬ 
nates (GBCs) (see |TU HU ESI 14S and references therein) which are defined below. 
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Fig. 2.1: Voronoi meshes of 8 x 8, 16 x 16, and 32 x 32 cells in Lloyd’s iteration. 



Fig. 2.2: History of Q a u, i and Q eq , i in Lloyd’s iteration for meshes with 8x8, 16 x 16 
and 32 x 32 cells. 


Definition 2.1. The generalized barycentric coordinates for a given n-gon T are 
the functions Xi : T— i = l,...,n satisfying 

(i) (Non-negativity) All Aj(x) for 1 < i < n are non-negative on T; 

(ii) (Linear precision) There hold 

n n 

YMx) = YM*)** : VxGT. 

i=l i= 1 


Take a set of GBCs on Tc and denote it by (Ai(^), A 2 (£),..., A n (^)). We can 













define a mapping Tt from Tc to T as 


x = JM£) = E A *(£) x - v ^ eT c- 

i—l 


Such a mapping is called a generalized barycentric mapping in literature and has 
important applications in several fields. An important non-trivial question is whether 
Tt defines a bijection or not. Fortunately, it has been answered positively for two 
types of GBCs on convex polygons, Wachspress coordinates nu and piecewise linear 
coordinates jT] [25]. 

The Wachspress coordinates [34} HE] are defined as 


Ai(x) 


Wi(x) . _ 

^ ' 7i ( \ ) ^ 1 , . . . , 


where 


Wi(x) = Xi_iXiX i+ i ]^[ XXjX J + i, 


i = 1,..., n 


and xyz denotes the signed area of the triangle with vertices x, y and z. The Wach¬ 
spress coordinates are linear for n = 3, bilinear for n = 4, and rational for general 
n. 

We now consider piecewise linear coordinates. Let Tt c be a triangulation of Tc- 
Let A i, i = 1,... ,n, be piecewise linear functions on Tt c satisfying A i(£j) = Sij and 
Aj = 1. This defines the piecewise linear GBCs associated with 7r c - Note that 
a different triangulation of Tc can lead to a different set of GBCs. Moreover, 7r c 
can be a triangulation using exactly the same vertices of Tq , or it can have extra 
vertices added as long as the extra vertices have their own and distinct barycentric 
coordinates specified. Furthermore, the piecewise linear bijection Tt from Tc to T can 
be obtained in a different but equivalent way: specify the same type of triangulations 
for Tc and T and then use piecewise linear mappings to map each individual triangle 
in Tc to the corresponding triangle in T. 

In practice, the piecewise linear barycentric mapping is much easier to compute 
than the Wachspress barycentric mapping. For this reason, we consider the former 
only in this work. In this case, the Jacobian matrix Jr of the mapping Tt is piecewise 
constant on each T £ T■ If we also use a piecewise constant approximation of M on 


T, from (2.9) and (2.10) we obtain the second set of mesh quality measures as 


_ trace((JT|jc) f MgJT|i<:) 

WaH ’2 “ T#T K&Tt 2det((J T |if) t M R -J T |^) 1 /2’ 


Q pn o = max max 
^ q ' T&TK&Tt 


det((Jr|ic) t Mr'J T | i< -)-\/det(Mrr) 


°7i,2 


where $t\k is the restriction of It on K , 
ah2 = iV^ 

trl TgTKgTt 

and N tr i is the total number of the triangles in the mesh. 


E E det((Jr|ic) t M/ 1 -Jr|^-)-\/ det(M^) 


(2.16) 

(2.17) 

(2.18) 
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As for Approximation 1, we test the second set of the mesh quality measures on 
the same Lloyd iterations shown in Fig. |2.1| Again, each n-gon is compared with a 
regular reference n-gon. There are two easy ways to define the triangular subdivision 
on T (as well as on Tc). (a) connect xi with all other vertices of T and get a 
triangulation; (b) connect all vertices to the arithmetic mean of vertices. Since T is 
convex, the arithmetic mean of its vertices lies inside T. We test both subdivisions. 
The results are reported in Figs. |2.3|and|2.4| 




Fig. 2.3: History of Q a u,2 and Q e q ,2 in Lloyd’s iteration for meshes with 8x8, 16 x 16 
and 32 x 32 cells. Here, the triangular subdivision (a) is used. 


Lloyds iteration for 8x8 mesh 





Fig. 2.4: History of Q a u,2 and Q eq ,2 in Lloyd’s iteration for meshes with 8 x 8 , 16 x 16 
and 32 x 32 cells. Here, the triangular subdivision (b) is used. 


Qaii,2 and Q eq ,2 have similar overall decreases, more dramatic drops at the first 
few iterations, and not monotone decreases as the mesh quality measures in Approxi¬ 
mation 1. They also provide correct indications of the behavior of Lloyd’s algorithm. 

On the other hand, one may notice that the values of Q a u,2 and Q eq , 2 are generally 
greater than Q a u,i and Q eq . 1 and more sensitive to quality changes (with more and 
stronger oscillations). This means that those Voronoi meshes have worse quality in the 
former’s eyes than in the latter’s eyes. The reason why Q a u,2 and Q eq ,2 are pickier 
is due to their construction: they favor triangulations that look good for the used 
piecewise generalized barycentric mapping. To explain this, we examine a random 
n-gon T in a CVT, as shown in Fig. |2.5| Although T is considered as of good shape 
by Lloyd’s algorithm, i.e., its generator is close to its barycenter, T may still contain 
short edges. When both the regular reference n-gon and T are cut into triangles which 
are compared one-by-one, the triangles associated with the short edges have a bad 
shape. 

It is also interesting to point out that subdivision (b) gives smaller Q a u,2 and 
Q eq ,2 than subdivision (a). 
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Fig. 2.5: Piecewise linear barycentric mapping from a regular reference hexahedron to 
T using triangular subdivision (b). Because of the short edges in T, some sub-triangles 
of T can have bad shape, which affects the value of Q a u, 2 - 


2.4. Approximation 3: employ a special design that uses infinitely 
many reference n-gons to ensure Tt be affine. This approximation is different 
from the previous two. Here, we introduce a set of infinitely many reference polygons 
of relatively good shape. Such a large set of reference polygons allows us to use affine 
mappings to define Tt- An important tool for the development is the connection 
between an arbitrary n -gon and a high dimensional simplex. Such a connection has 
been known in geometry and discrete mathematics (e.g., see [49l Chapter 0]). Here 
we use it to analyze polygonal mesh quality. 

Consider a set of GBCs A,;, i = 1, ...,n, on T G T- Recall that 


x = ^AiXi, VxeT (2-19) 

i= 1 

where A = [Ai, A 2 , • • • , A ra ] 4 lies in the set 

n 

S n -1 = {[Ai, A 2 ,..., A nY y> = l and A* > 0, 1 < i < n}. 


The set S n _ 1 is a regular (n— l)-simplex in K", meaning that it has complete symmetry 
with regard to transitions over vertices, edges, and higher dimensional faces, etc. S n -± 
is contained in an (n — 1)-dimensional hyperplane ]C" =1 Ai = 1 in R n and is essentially 
an (n — 1)-dimensional simplex living in the n-dimensional space. 


Equation (2.19) can be viewed as a linear mapping from S n -1 onto T. We shall 


add a translation to this mapping to make the center of mapped to the origin, 
which will greatly simplify the derivation. Denote the arithmetic mean of all vertices 
of T, which will be called as the arithmetic center of T, by 


1 " 

xt = y>, 

n L ^ 


2 — 1 


The arithmetic center of a convex polygon always locates inside the polygon. Equation 


(2.19) can be rewritten into 


x - x T = ^2 Aj(x, - x T ) 


( 2 . 20 ) 


Denote by x = x - xt a translation in coordinates which moves xt to the origin. 
Then, (2.20) can be viewed as a linear mapping from A to x, with the geometric 
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center of S n - 1 , denoted by A c = ^] 4 , being mapped into the origin x = 0. 

For simplicity and with abusing of the notation, it is convenient to still denote the 
coordinate for T by x, which should actually be understood as x, together with the 
assumption that = 0. 

The mapping from S n -i to T can be written as 


S n _i 


T, 


’Ft(A) = Bt^ — [x x X 2 


Clearly, one has 


0 

0 


Xt = 4't(A c ). 


X„] A. 


( 2 . 21 ) 


When n = 3, (2.21) together with the constraint 
solvable 3x3 linear system for non-degenerate T, 


A» = 1 gives a uniquely 


Bt 

l 4 


A4 


Xx 

1 


x 2 X 3 
1 1 


A = 


i.e., the mapping H/t is invertible. But for n > 3, the mapping Tt is not invertible 
as there is a non-trivial kernel. This can also be explained by the fact that T is a 
2-manifold while S n -\ is an (n — l)-manifold. 

To understand the linear mapping 'Ft from S n -\ to T, we consider the singular 
value decomposition of the matrix B t € R 2xra , 


F?t — C/t^tIt: 


where 


Ut 


Vt 


[Ui ; T U-2,t] G M 2x2 , 


[ V 1,T V 2,T • ' ' v «,t] € 



0 0 
(72,T 0 


p2xn 


<7\,t > <72,t are the singular values, and Ut and Vt are orthogonal matrices. The 
singular values are non-zero as long as the polygon T does not degenerate into a line 
segment. We further decompose 


(7l,T 

0 


0 

(72,T 



0 

1 


0 ••• 0 

0 ••• 0 


= D t Qt- 


Then, the singular value decomposition of Bt can be rewritten into 


Bt ~ (UtD t Ut){UtQtVt) — Jt Pt- 


Now, the geometric meaning of the linear mapping Tt (with coefficient matrix Bt) 
can be explained as follows (see an illustration given in Fig. 2.6). 

First, the linear mapping Tt applies to S n -1 with a linear transformation speci¬ 
fied by 


Pt = UtQtVt = Ut 
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Since QtVt is an orthogonal projection from S n -i onto the plane spanned by v^t 
and v 2j t, Pt is an orthogonal projection followed by a 2D rotation/reflection Ut- 
The result is an n-gon 


At — Uj 


'1 ,T 

I 

V 2,T 


(Sn-l)- 


( 2 . 22 ) 


Second, the mapping applies to K T with a linear transformation Pt($,) = Jt£ = 
UtDtU resulting in T = Ft{Kt )■ Note that It is symmetric and positive definite. 
Thus J~t is an anisotropic scaling , i.e., scaling by factors a± t T and 02,t in the directions 
of Ui t and U 2 ,Tj respectively. Polygons T and Kt are affine similar under the 
anisotropic scaling Ft- 



Fig. 2.6: Illustration of the linear mapping from simplex S n -\ to a polygon, which is 
decomposed into two affine mappings. 


In the following we show that Kt has several nice properties. First of all, because 
T and Kt are affine similar to each other, we have the following three propositions. 

Proposition 2.2. T is convex if and only if Kt is convex. 

Proposition 2.3. T is non-degenerate if and only if Kt is non-degenerate. 
Proposition 2.4. The arithmetic center of Kt locates at the origin. 


More importantly, Kt has a relatively good shape and can serve as a reference 
n-gon. To show this, we first need to clarify what it means for an n-gon to have a 
relatively good shape. For triangles, it is common to use the ratio between the radii 
of its inscribed circle and circumscribed circle as a measure. But an n-gon with n > 4 
may not have inscribed or circumscribed circles. Therefore, we introduce the following 
definition. 


Definition 2.5. For a given polygon T, its in-radius is the maximum radius of 
all disks contained inside T , and its outer-radius is the minimum radius of all disks 
containing T. 


Proposition 2.6. LetT be a convex n-gon and Kt be defined in (2.22). Then, 
the in-radius of Kt is greater than or equal to an d the outer-radius of Kt 


is less than or equal to y Consequently, the ratio p between the in-radius and 

the outer-radius of Kt is bounded by 


1 


n — 1 


<P< 1, 


which depends on n but not on the shape of Kt- 
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Proof. Recall that ,S' n _i lies inside the hyperplane A?: = 1- By Proposition 
|2.4| and linearity, one has 


0 

0 


TV A,. = Hr 


Vl,T • A c 

V 2,T ’ A c 


which implies that • A c = V 2 ,t ■ A c = 0. Consequently, vectors in span-jV^T, V 2 ,t} 
are orthogonal to A c = [A,..., A] 4 , which is also the normal direction of the hyperplane 
Sr=i = 0. In other words, the 2D plane spanned by v 1>T and v 2 ,t lies inside the 
hyperplane = which is parallel to the hyperplane XliLi A* = 1- This 

is important as it guarantees that the orthogonal projection of any ball inside the 
hyperplanc J]”=i A» = 1 into spanjvijy V 2 ,t} must be a circular disk, as illustrated 
in Fig. [2~7l 


Hyperplane 



Projection plane 



Projection plane 


Fig. 2.7: Orthogonal projection from a hyperplane to a 2D plane. Left: If the 2D 
plane lies inside a parallel hyperplane, then projection of any ball in the hyperplane 
gives a circular disk; Right: if not parallel, the projection becomes an ellipse. 


By Proposition |2.2[ polygon Kt is also convex. Hence the projection Pt maps 
the inscribed ball of S n -\ into an inner disk of Kt, and the circumscribed ball of S n —i 
into an outer disk of Kt- One simply needs to compute the radii of the inscribed and 

the circumscribed ball of S n _ i, which are just ^/ n ( n 1 _ 1 ^ and This completes 

the proof of the lemma. □ 


Remark 2.7. For n > 3, the projection of the inscribed (circumscribed) ball of 
S n — i under Pt is not necessarily the largest disk in Kt (the smallest disk containing 
Kt)- 

Because of the uniform bound for p stated in Proposition |2.6[ K T has a relatively 
good shape and can thus be used as a reference n-gon. 


A few possible AV’s are shown in Fig. |2.8[ When n = 3, the simplex S 2 is an 
equilateral triangle lying in the plane Ai + A 2 + A 3 = 1 and the condition vi • 1 = 
v 2 • 1 = 0 implies that span{vi,v 2 } is just the plane Ai + A 2 + A 3 = 0. Thus the 
reference 3-gons, or the reference triangles, ar e sim ply rotated/reflected equilateral 


triangles with edge length \/2, as shown in Fig. 2.8 In other words, the (equilateral) 


reference triangle is unique up to rotation/reflection. When n > 3, on the other hand, 
the reference n-gons will no longer be unique even after rotation/reflection. However, 
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Fig. 2.8: Examples of reference polygons. The dotted circles are the projection of the 
inscribed and circumscribed balls of S n -1 under Pt- 


they should all lie between two circles with radii y n(n-i) an d y as stated in 
Proposition |2.6[ and they should all be convex and non-degenerate, as verified in 
propositions |2.2| and |2.3| 

Now we are ready to define the mesh quality measures using the reference n- gons. 
Recall that each T has a unique Kt and Tc is the collection of all those Kps. (We do 
not need to compute Kt 's when evaluating the mesh quality measures.) The Jacobian 
matrix It of Pt is defined as 


Jy = Ut 


01,T 
0 


0 

02, T 


Uf 


T) 


(2.23) 


where Ut and <ti,t and <T 2 ,t are the left singular vectors and singular values of the 
matrix Bt = [xi — xr,...,x n — xy], (Here we put back the translation just for 
convenience.) Using the metric tensor approximation (2.1), we can then define the 
mesh quality measures as 


where ly is defined in (2.23), 


N, 


trace(JyMrl7’) 
tJt 2det(jyM T J T ) 1 /2’ 

(2.24) 

det(iyMTlT) -J det(M-r) 
max-, 

(2.25) 

Ter 3 


— 'y ' det(I^M tIt) det(M t), 

(2.26) 


Ter 


and N p is the number of the polygons in T. 

We have tested the third set of mesh quality measures numerically on the same 
Lloyd iterations shown in Fig. 2.1 The results are reported in Fig. |2.9| Comparing 
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with the numerical results for Q a u,i and Q eq ,i (cf. Fig. 2.2 and Table 2 . 11 , one may 


notice that Q a u,3 is slightly smaller than Q a u,i but Q eQ} 3 is slightly larger than Q eq 1 
for the same mesh. This is because in the computation for the first approximation, 
regular n -gons are used as reference n-gons, which gives a precise control of the size of 
reference n- gons; while in the third approximation, the size of reference n-gons varies 
a bit as can be seen in Fig. 2.8 Overall, Q a u ,3 and Q eq ,3 also provide a good measure 


on the quality of Voronoi meshes. A major advantage of this approximation is that 
all TV’s are affine, which will be useful in error estimation in numerical solution of 
partial differential equations. 



Number of iterations Number of iterations Number of iterations 


Fig. 2.9: History of Q a u, 3 and Q eq ,3 in Lloyd’s iteration for meshes with 8 x 8 , 16 x 16 
and 32 x 32 cells. 


2.5. Summary. The three sets of mesh quality measures discussed in the previ¬ 
ous subsections can more or less be viewed as discretizations or numerical approxima¬ 
tions of the continuous level alignment and equidistribution quality measures defined 
in (2.10). Moreover, they all adopt the idea of evaluating the quality of a polygo¬ 
nal mesh by comparing its elements to their counterparts in the computational mesh 
which can be a conventional mesh or a collection of reference polygons. Numerical 
experiment shows that they all provide correct measures for the quality of Voronoi 
meshes generated by Lloyd’s algorithm. Some of the main features of these mesh 
quality measures are elaborated in the following. 

• Qau, 3 (2.241 and Q eq , 3 ( |2.25 ) are fully determined by the mesh T. Specifically, 
the computational mesh 7c is a collection of Kt’s which are shown to have 
reasonably good quality and determined by the coordinates of the vertices 
of T. The mapping Tt is an affine mapping from I\t to T, which is also 
determined by the coordinates of the vertices of T. 

• Qau, 1 (2.13) and Q eq , 1 (2.14| are not fully determined by the mesh T. The 
computational mesh 7c can be chosen by the user to be a mesh or a collection 
of reference n-gons. The measurement of the quality of T depends on the 
choice of 7c- The mapping Tt is approximated by an affine mapping that is 
obtained by least squares fitting to the correspondence between the vertices 
of T £ T and Tq £ 7c and thus fully determined by T and 7c- 

• Qau, 2 (2.16) and Q eq ,2 (2.17) are not fully determined by the mesh T. The 
computational mesh 7c can be chosen by the user to be a mesh or a collection 
of reference n-gons. The mapping Tt is specified using generalized barycen- 
tric mappings such as Wachspress and piecewise barycentric coordinates. The 
measurement of the quality of T depends on the choice of 7c as well as the 
choice of Tt- The tests using Voronoi meshes generated by Lloyd’s algorithm 
(cf. Table 2 . 1 ) suggest that Q a u,2 and Q eq ,2 provide the toughest measures 
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(which give biggest values and are hardest to satisfy) among the three sets of 
mesh quality measures. 


Table 2.1: Q a u and Q eq at selected Lloyd’s iteration steps for the 32 x 32 mesh. 
Subdivision (b) was used for Q a u,2 and Q eq ,2- 


Iter. 

Qali, 1 

Qeq, 1 

Qali, 2 

Qeq,2 

Qali,3 

Qeq,3 

0 

3.4362 

3.0130 

3.4362 

7.3202 

3.4362 

3.6641 

2 

1.8519 

2.5239 

2.9940 

3.3286 

1.8851 

2.6679 

8 

1.3927 

1.5355 

2.9190 

2.5647 

1.3869 

1.9022 

43 

1.1394 

1.3771 

2.7847 

2.1333 

1.1370 

1.4703 


3. Anisotropic polygonal mesh adaptation. In this section, we study the 
anisotropic adaptation for polygonal meshes through a moving mesh method based on 
the so-called MMPDE (moving mesh PDE) [25] (25]. Notice that a mesh is completely 
determined by two data structures: the coordinates and connectivity of vertices that 
form polygons. In a certain range, one can move the vertices without changing mesh 
topology or tangling the mesh. The moving polygonal mesh method uses this idea 
and implements it in an iterative manner. 

1. Initialization: Given an initial physical mesh T < - 0 ' 1 for S2; 

2. Outer iteration (fc = 0,1,...): 

(a) Update the metric tensor Mb fc ) based on the information available at 
the current iteration. The information includes the current mesh 
and the physical solution that is obtained by solving the underlying 
PDE on the current mesh T^; 

(b) Find a way to move the vertices of the physical mesh so that the new 
mesh T (fe+1) has a better quality under the metric 

We now give a detailed description of Step[2b] The moving mesh strategy is closely 
related to Q a u,2 and Q eq , 2 (with piecewise barycentric mappings and subdivision (b)). 
Recall that they are the toughest set of mesh quality measures among the three, and 
it is reasonable to expect that the others can be satisfied when they are satisfied. 
Moreover, Q a u,2 and Q eq ,2 are associated with a triangulation of each polygon in 
the mesh. One advantage of using piecewise barycentric mappings to define local 
mappings Tt is that it allows us to reuse part of the code previously developed in 
[24| for moving triangular meshes. The triangular subdivisions of all T £ T used 
in defining Tt, combined together, can be viewed as a triangular mesh itself. On 
the other hand, we should emphasize that the polygonal moving mesh algorithm is 
substantially different from a triangular moving mesh algorithm since the former aims 
to optimize each polygon whereas the latter focuses on the quality of each triangle. As 
shown in Fig. |2.5[ some of the sub-triangles of a good polygon can indeed have a bad 
shape. 

We assume that a reference computational mesh Tc has been chosen for the mesh 
movement purpose. For example, we can use a CVT generated by Lloyd’s algorithm. 
We denote the triangular meshes resulting from the triangulation associated with 
Qau,2 and Q eq ,2 for Tc and 7 ^ fc+1 ) by Tf c and 77-1*+i>, respectively. Consider the 
function 

4({£j,{xf +1) })= £ \K\G(J K ,det(J K )M K ), (3.1) 

ifeT T( fc + i) 
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where {£,} and {x( fc+1) } denote the coordinates of the vertices of Tf c and Vf(k+ p, 
respectively, M k is the average of on K, Jk is the inverse of the Jacobian matrix 
lie of the affine mapping from Kc £ Tf- to K £ Tq-ik+i) , and 


G(J k , det( J k ),M k ) = ^ \/ det ( M ^) (trace( J k ^&k j k)Y 


] V det(Mjf) ( d f ( ^M 


V / det(M if )) 


(3.2) 


The function (3.1) is a discretization of a continuous meshing functional EDI 122 and 
it is known that minimizing Ih will tend to make the mesh to satisfy the alignment 
and equidistribution conditions associated with M( fc ) and thus to have a better quality 
under the metric 

Since Tf c is known, the new mesh {x| l+1 *} or Tj-ik+i) can be obtained either by 

directly minimizing Ih. or by first differentiating Ih with respect to {x,- ,:+l } and then 
solving the corresponding algebraic equations. In either case, an iterative method has 
to be used because Ih is highly nonlinear. The iterative solution of the new mesh 
forms the inner iteration. 

Recall that is defined on ■ During the inner iteration, the metric tensor 
needs to be updated constantly (through interpolation) on approximate meshes of 
T (fc+1) that will have different vertex locations than . This can be expensive even 
if linear interpolation is used. 

To avoid this difficulty, an indirect method for finding 77 - 0 + 1 ) is used. We replace 
(3.1) with {x(^}) (the same form but with difference meshes), where {£,} 


denotes the vertices of a computational mesh 7j- c and {x- i '} denotes the vertices of 
the physical mesh 7)-<&>. We then minimize the new function with respect to {£;} while 
{x( fe -*} remains fixed. An advantage of this is that, since is fixed during the inner 
iteration, there is no need for updating Moreover, minimizing 7/i({£j}, {x-^}) 

also tends to make the mesh to satisfy the alignment and equidistribution conditions 
associated with Thus, it is reasonable to expect that the correspondence 


{xi*>} = $„({&}) 


{xj fc+1) } = *„({&}), 


is close to the correspondence 


which is determined by minimizing ( |3.1[ ). Then, we can define the new mesh as 

{xf +1) } = <!>,({&}), 


which can be computed by linear interpolation as is completely determined once 
{x- fe ^} and {£,} are known. 

Now, the inner iteration is viewed as to find the steady state of a modified gradient 
flow of the function Ih({€i}, {x.- fc1 }). The mesh equation is recorded below without 
derivation. The interested reader is referred to [22 for detailed derivation. 


= ^ E \K\v. 


K 


Keui 


k €i(0) = &, 
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i = 1,..., N v , f > 0 
i = 1,... ,N V 


(3.3) 








where N v is the number of vertices of Tj-(k) , LOi is the patch of triangles associated with 
vertex xin T-ytk) , %k is the local index of x in K, vf K is the local mesh velocity 
associated with the vertex of K, r > 0 is a constant parameter used to adjust the 
time scale of mesh movement, and P = (Pi,... ,Pn„) is a positive function used to 
make the mesh equation to have desired invariance properties. The local velocities 
are given by 


( 0*1 


= -E 


K 


dG 

DJk 


dG det (E Kc ) , 
ddet (Jk) det(-Eff) Kc ’ 


,.k 


i=i 


,K 


(3.4) 


where 

^ \/ det(Mif) trace( J K Mj} j£) M^J^, 

dG 8 det (Jk) 
ddet(J K ) 3 y/det(M K )' 


The balancing function in (3.3) is chosen to be Pi = det(M(x-^))5 such that (3.3) is 
invariant under the scaling transformation M —> cM. 

The mesh equation (3.3) should be modified properly for boundary vertices. For 
fixed boundary vertices, the corresponding equations should be replaced by 


dti 

dt 


= 0 . 


For boundary vertices on a curve represented by <fi(Z) = 0, the corresponding mesh 
equations should be modified such that its normal component along the curve is zero, 
i.e., 




dZi 

dt 


= 0. 


The mesh equation (3.3), along with proper modification for boundary vertices, 


can be solved by any ODE solver. (Matlab’s ODE solver odel5s is used in our com¬ 
putation.) In principle, it should be integrated until a steady state is obtained. Since 
finding just represents one step of the outer iteration, we integrate (3.3) only up 
to t = 1 (with r = 1/300) to save CPU time. 

We now discuss the choice and computation of the metric tensor. We choose the 
metric tensor based on optimizing the L 2 norm of error for piecewise linear interpola¬ 
tion [221 [2Z3 - The main reason for this choice is that it is simple, problem independent, 
and effective. It reads as 


= det (a h I + \H(u h )\) 6 [a h I + \H(u h )\\ , 


(3.5) 


where Uh is an approximate solution, H(uh ) is the recovered Hessian of Uh, \H(uh)\ 
is the eigen-decomposition of H(uh ) with the eigenvalues being replaced by their 
absolute values, and the regularization parameter ah > 0 is chosen such that 

I \J det(M) dx = 2 f det (\H(uh)\) ^dx. (3.6) 

Jn J n 

It has been shown in [2Sj that when the recovered Hessian satisfies a closeness as¬ 
sumption, a linear finite element solution of an elliptic boundary value problem on 
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a simplicial mesh computed using the moving mesh algorithm converges at a second 
order rate as the mesh is refined. Numerical examples presented later show that the 
same strategy seems to work well for polygonal meshes too. We use a Hessian recovery 
method based on a least squares fit. More specifically, a quadratic polynomial is con¬ 
structed locally for each vertex via least squares fitting to neighboring nodal function 
values and an approximate Hessian at the vertex is then obtained by differentiating 
the polynomial. 

Finally, we point out that currently there is no built-in mechanism in the al¬ 
gorithm to prevent mesh elements from becoming non-convex, or to prevent mesh- 
tangling. Such issues, as well as the convergence of the algorithm, are generally 
difficult to analyze theoretically |5|. On the other hand, the continuous meshing 
functional corresponding to Ih({£i}> {x| fc) }) is coercive and polyconvex and has a 
minimizer [2U1 . Moreover, numerical results to be given in Section [4] show that the 
algorithm is efficient and robust. 

4. Numerical results. In this section, we present numerical results obtained 
with the moving mesh algorithm given in Section [3] for Laplace’s equation -A« = / 
in = (0, l) 2 subject to the Dirichlet boundary condition. Two examples with the 
following exact solutions are tested, 

Example 1: u = tanh(40y — 80x 2 ) — tanh(4Cte — 80 y 2 )\ 

Example 2: u= y/0.5 (r — x) — 0.25r 2 , r = \J x 2 + y 2 . 

These examples are solved on polygonal meshes using the Wachspress finite element 
method [46] . an H 1 conforming finite element method using the Wachspress barycen- 
tric coordinates as basis functions. It is known that for smooth exact solutions and 
sufficently fine shape-regular quasi-uniform polygonal meshes, the Wachspress finite 
element method has the asymptotic convergence order 0(h) in H 1 semi-norm and 
0(h 2 ) in L 2 norm, where h is the characteristic size of the mesh. If one considers a 
quasi-uniform polygonal mesh with N x N polygons, the characteristic mesh size h. 
roughly equals iV -1 . Consequently, the optimal asymptotic order of the approxima¬ 
tion error for the Wachspress finite element method can be expressed into 0(N~ 1 ) in 
H 1 semi-norm and 0(N~ 2 ) in L 2 norm. 



0 0.5 1 

Fig. 4.1: Contour plot of the exact solution for Example 1. 


Let us briefly explain the reason why we pick these two examples. For Example 1, 
the value of u ranges from —2 to 2 in Q. An equidistant contour plot of u is shown in 
Fig. 4.1 Apparently, u changes rapidly near the curves y — 2x 2 = 0 and x — 2 y 2 — 0, 


20 









while it is almost constant everywhere else. It exhibits a strong anisotropic behavior 
in the gradient and Hessian near those curves. When a uniform isotropic mesh is used 
to discretize Example 1, the number of the elements has to be very large to resolve 
the rapid changes of u near the curves. On the other hand, less elements can be 
used for the same accuracy when using an adaptive anisotropic mesh, which is dense 
around the regions with the rapid changes of u and coarse in places where u is almost 
constant. Here, we will show that the moving mesh algorithm described in Section [3] 
is capable of generating high quality anisotropic adaptive polygonal meshes optimized 
for this problem in terms of the L 2 norm of the approximation error. 

Example 2 is a well-known example with a corner singularity at (0,0), and the 
exact solution u is in H^~ e {VL) for arbitrarily small e > 0. When discretized using 
quasi-uniform meshes, the best approximation error that can be reached has asymp¬ 
totic order O(h 0 - 5 ) = 0(N ~°- 5 ) in H 1 semi-norm and 0(/i 15 ) = 0(N~ 1 - 5 ) in L 2 norm. 
It is emphasized that, no matter how fine the uniform mesh is, the asymptotic order 
of approximation error cannot be improved due to the intrinsic low regularity of the 
exact solution. Later on we will show that adaptive meshes generated by the moving 
mesh method not only can lead to smaller errors but also improves the convergence 
order to the optimal 0(N~ 2 ) for the L 2 norm. 

For both examples, we set flc = fl and use CVTs generated on £l c by Lloyd’s 
algorithm as the reference computational mesh Tc- According to the numerical results 
given in Section [ 2 ] the mesh Tc has good quality under the Euclidean metric. We 
take the initial physical mesh T® = Tc- 

We start from a reference computational mesh Tc with 32 x 32 polygons and 
apply the moving mesh algorithm with 10 outer iterations to Example 1. The initial 
mesh and the physical mes hes after 1 and 10 outer iterations, i.e., T*- 0 ), T*- 1 -* and 
T (10) , are shown in Fig. 


4.2 


O ne ca n see that the meshes correctly capture the rapid 
changes of the solution (cf. Fig. 4.1). In addition, a close view of T*- 10 -* near (0.5, 0.5) 
(shown in Fig. 4.2) clearly shows the anisotropic behavior of the mesh elements. The 
history of alignment and equidistribution measures is reported in Fig. |4.3[ Here, 
Qaii, 1 , Qaii, 2 and Qeq, l, Qeq, 2 are computed by comparing the physical mesh with the 
reference computational mesh Tc, while Q a u, 3 and Q e q, 3 are computed by comparing 
each polygon T in the physical mesh with its own reference Kt- One can see that the 
moving mesh algorithm reduces the mesh quality measures although the reduction is 
not monotone. Moreover, Q a u, 2 is slightly bigger than Q a n,i and Q a u, 3 and Q eq , 3 
is slightly bigger than Q eq ,i■ These are consistent with what we have observed in 
Section [2] Fig. |4.3| also shows that all three sets of mesh quality measures have very 
similar evolution patterns. 



Fig. 4.2: Example 1, mesh with 32 x 32 cells. From left to right: T^°\ T^\ T^ w \ a 
close view of Tnear (0.5, 0.5). 
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Fig. 4.3: Example 1, mesh with 32 x 32 cells. History of Q a u and Q eq . Subdivision 
(b) was used for computing Q a u,2 and Q eq ,2- 


L2 norm of the error 



HI semi-norm of the error 



Fig. 4.4: Example 1, a mesh with 32 x 32 cells. History of L 2 norm and H 1 semi-norm 
of the error u — Uh is plotted as a function of the outer iteration number. 


A more important question is whether these adaptive meshes actually help reduce 
the approximation error or not. To examine this, we computed the L 2 norm and the 
H 1 semi-norm of u — u^, where Uh is the finite element solution using the Wachspress 
finite element method, on these adaptive meshes. The results are reported in Fig. |4.4| 
The figure clearly shows the effectiveness of the moving mesh method in reducing the 
approximation error, as both the L 2 norm and the H 1 semi-norm of the error are 
reduced for about 10 times after applying 10 outer iterations. Interestingly, although 
the algorithm does not reduce Q a u or Q eq monotonically, it seems to reduce the L 2 
norm and the H 1 semi-norm of the error monotonically for Example 1. We also notice 
that the majority of reduction occurs within the first few outer iterations. 

With the above observations, we continue testing the method for Example 1, on 
meshes with different numbers of cells. Consider meshes with N x N polygonal cells, 
for N = 8, 16, 32, 64, and 128. Again, the reference computational meshes are taken 
as CVTs generated by Lloyd’s algorithm. Optimal rates of convergence, 0(N^ 1 ) in 
H 1 semi-norm and 0(N~ 2 ) in L 2 norm, usually cannot be achieved when the mesh 
is not fine enough to resolve all detail of the solution. In Table |4.1| it can be seen 
that on quasi-uniform mesh asymptotic order of errors in H 1 semi-norm is less 
than 0(N~ 1 ) and improves as N increases. The asymptotic order of the error in L 2 
norm on T*- 0 - 1 appears to be larger than 0(N~ 2 ) at the beginning, which is indeed an 
indication that very coarse meshes cannot fully resolve the detail of the solution and 
thus result in bad approximations. Again, when N increases, the asymptotic order of 
the error in L 2 norm improves until it reaches 0(N~ 2 ). 
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L 2 norm of the error H 1 semi-norm of the error 




Fig. 4.5: Example 1, mesh with TV x TV cells, for TV = 8, 16, 32, 64, and 128. The L 2 
norm and the H 1 semi-norm of the error u — Uh after different numbers of MMPDE 
outer iterations are plotted as functions of TV. 


Table 4.1: Example 1, mesh with TV x TV cells, for TV = 8, 16, 32, 64, and 128. The 
L 2 norm and the H 1 semi-norm of the error on T-°\ i.e., no MMPDE iteration, and 
T^ 5 \ i.e., 5 MMPDE iterations. Here, the asymptotic order of the error is computed 
using two consecutive meshes with respect to L ■ 



On mesh 

On mesh T (5) 


L 2 norm 

H 1 semi-norm 

L 2 norm 

H 1 semi-norm 

TV 

error 

order 

error 

order 

error 

order 

error 

order 

8 

1.50e+0 


1.63e+l 


4.15e-l 


1.16e+l 


16 

2.16e-l 

2.8 

1.13e+l 

0.5 

2.65e-2 

4.0 

4.17e+0 

1.5 

32 

5.45e-2 

2.0 

7.32e+0 

0.6 

3.54e-3 

2.9 

1.51e+0 

1.5 

64 

1.63e-2 

1.7 

4.25e+0 

0.8 

1.Ole-3 

1.8 

7.04e-l 

1.1 

128 

4.39e-3 

1.9 

2.23e+0 

0.9 

2.73e-4 

1.9 

3.51e-l 

1.0 


Both the H 1 semi-norm and L 2 norm of the approximation error are significantly 
improved by the anisotropic adaptive mesh algorithm. In Fig. |4.5| the approximation 
error is reported for different mesh sizes and different numbers of MMPDE outer 
iterations. We can see that the approximation error can be effectively reduced by 
performing just a few MMPDE outer iterations. To better compare the quantity of 
the error, we also list the results o n th e initial mesh and on the physical mesh after 5 
MMPDE outer iterations in Table 


4.1 


Asymptotic orders in terms of A are reported 
in the table too. Improvements in both the H 1 semi-norm and the L 2 norm and in 
convergence order can be observed. 

Next, we test Example 2 under the same settings as for Example 1. We start 
from an initial mesh with 32 x 32 cells and apply the MMPDE algorithm. The initial 
mesh and physical meshes after 1 and 10 outer iterations of MMPDE are shown 

and the L 2 norm and H 1 semi-norm of 
One may notice that Q a u, i, 


in Fig. 


4.6 


The history of Q a u and Q, 


the error is reported in Figs. 


4.7 and 4.10 


eq 


and 


Qau, 3 stay almost constant while Q eq , l, Qeq, 2 , and Q eq , 3 increase slightly during the 
MMPDE iterations. We first point out that the corner singularity in this example 
is essentially isotropic. Therefore, any isotropic mesh has good alignment under any 
metric based on the recovered Hessian of the solution. As a result, we expect that Q a u 
remains small and constant during the MMPDE outer iterations. To better illustrate 
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this, in Fig. 4.8 we plot the piecewise constant functions q a u,i 
Approximation 1, i.e., 


and 


Qeq , 1 


defined in 


Qali,l — 


Qeq ,1 — 


tr ace ( [Et Etp c (Et c E^ c ) “ 1 ] 4 Mt [Et Etp c {Et c Etp c ) -1 ]) 
2det({E T E t Tc (E Tc E t Tc )^}m T {E T E t Tc (E Tc E t Tc )-^ , 
det ( Et E^ c ( Et c Etp c ) “ 1 )\J det(My) 


where o> tl is defined as in (2.15). The ranges of q a u,i and q eq ,\ are [l,oo) and (0, oo), 
respectively. Notice that Q a ii,i and Q eq .i are the maximum norm of q a u,i and q eq ,i- 


The left panel of Fig. 4.8 shows that q a u, i stays almost constant on the entire domain 


and its distribution is independent of the corner singularity. (We plot only q a u t i and 
q eq ,i in Fig. 4.8 since the corresponding quantities for other mesh quality measures 


behave similarly.) 

On the other hand, we recall that the exact Hessian for this example is infinite at 
(0, 0). Thus, the more adaptive the mesh is, more elements are moved toward the ori¬ 
gin and the computed Hessian has greater values there and overall it is becoming less 
smoother. As a consequence, it is harder to generate a mesh with perfect equidistribu- 
tion and we see that Q eq is getting bigger as the MMPDE outer iterations. Moreover, 
Q eq is defined as the maximum norm of q eq , which tells about the mesh quality of 
worst polygons. As shown in the right panel of Fig. |4.8[ worst polygons occur near 
the origin. The L 2 norm of the quantities is shown in Fig. 4.9 The results show that 
the L 2 norm decreases almost monotonically, indicating that the MMPDE algorithm 
has successfully improved the quality of the majority of mesh elements although the 
worst polygons remain “bad”. 


The L 2 norm and H 1 semi-norm of the approximation error are shown in Fig. 4.10 


Both of them decrease monotonically for Example 2, which further suggests that the 
MMPDE algorithm works effectively for Example 2. 



Fig. 4.6: Example 2, mesh with 32 x 32 cells. From left to right: T^°\ T^ w \ a 
close view of near the origin. 


Finally, we test the MMPDE algorithm for Example 2 on N x N meshes with 
different N. As mentioned earlier, the error u — u/, has at best the asymptotic order 
0(N~ 0 ' 5 ) in H 1 semi-norm and 0(N~ 1 - 5 ) in L 2 norm on Nx N quasi-uniform meshes. 
This can be seen clearly in Table fO| In Fig. 4.11 and Table 4.2 the error for Example 


2 under different mesh sizes and MMPDE outer iterations is reported. From Fig. 4.11 


it is clear that increasing the number of MMPDE outer iterations affects the order 
of L 2 norm and H 1 semi-norm of the approximation error. The numerical values 
of the asymptotic order reported in Table |4.2| give a clearer comparison. After 5 
MMPDE outer iterations, the L 2 norm of the approximation error achieves almost 


24 































MMPDE outer iter. 



Fig. 4.7: Example 2, mesh with 32 x 32 cells. History of Q a u and Q eq . Subdivision 
(b) was used for Q a u, 2 and Q eq , 2 - 




Fig. 4.8: Example 2, distribution of q a n, 1 and q eq ,i on mesh with 32 x 32 cells after 
10 MMPDE outer iterations. 


0 (N~ 2 ). The asymptotic order of H 1 semi-norm also improves, but not to the optimal 
0 (N~ 1 ) order. This may be because the algorithm is designed to optimize the L 2 
norm of the error instead of the H 1 semi-norm of the error. To see this, we use the 
metric tensor 


= det ( a h I + \H(u h )\) 4 \\a h I + \H(u h )\p [a h I + \H(u h )\\, 


(4.1) 


where ah is determined through the equation ( |3.6| ). This metric tensor is based on 
minimizing the H 1 semi-norm of linear interpolation error [22] , The numerical results 
are shown in Fig. 4.12| and Table 4.3 The convergence order in the H 1 semi-norm 
improves (around 0.9) while maintaining the second order rate for the L 2 norm. 

5. Conclusions. In the previous sections we have studied anisotropic mesh qual¬ 
ity measures and adaptation for polygonal meshes. Three sets of alignment (for shape) 
and equidistribution (for size) quality measures have been developed. They can be 
viewed as numerical discretizations or approximations of continuous-level measures 
(2.10). Their major features are summarized in subsection 2.5 Numerical examples 


show that they all can give good indicators for the quality of polygonal meshes. 

Among the three sets of measures, Q a ii, 2 (2.16) and Q eq ,2 (2.17) provide the 


toughest measurement. One of the special cases for defining Q a u,2 and Q eq ,2 is to 
use piecewise linear generalized barycentric mappings which can be based on trian¬ 
gulation of each polygonal cell into triangles. The collection of those triangles forms 
a triangular mesh, through which many adaptive mesh algorithms designed for tri- 
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Fig. 4.9: Example 2, mesh with 32 x 32 cells. History of ||g a /i||L 2 and ||<7e 9 ||L 2 - Sub¬ 
division (b) was used for Q a u,2 and Q eq ,2- 



Fig. 4.10: Example 2, mesh with 32 x 32 
of the error u — u^. 


HI semi-norm of the error 



. History of L 2 norm and H 1 semi-norm 


angular meshes can be adopted for polygonal mesh adaptation. Though special care 
needs to be taken since a good polygonal cell can have bad triangles and algorithms 
should focus on improving the quality of polygonal cells instead of triangles. 

Along this line, an anisotropic adaptive polygonal mesh algorithm based on the 
MMPDE moving mesh method has been devised. For a given reference computational 
(polygonal) mesh, it moves the mesh vertices such that the adaptive polygonal mesh 
has a good quality with reference to the reference one. In the numerical examples, 
the reference computational mesh has been taken as a CVT generated with Lloyd’s 
algorithm. Numerical results confirm that the proposed anisotropic adaptive polygo¬ 
nal mesh algorithm is able to produce desired mesh concentration and lead to a more 
accurate solution than using a uniform polygonal mesh. 

Finally, we would like to point out that the main idea in this work can be gen¬ 
eralized to three-dimensions to define anisotropic polyhedral mesh quality measures 
and generate anisotropic adaptive polyhedral meshes. However, the technical details 
of the implementation and the effectiveness and robustness of the method in three- 
dimensions remain to be examined in future work. 
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Vr norm of the error 


H 1 semi-norm of the error 




Fig. 4.11: Example 2, mesh with TV x TV cells, for TV = 8,16,32,64,128. The L 2 
norm and H 1 semi-norm of the error u — Uh after different numbers of MMPDE outer 
iterations are plotted as functions of TV. 


Table 4.2: Example 2, mesh with TV x TV cells, for TV = 8, 16, 32, 64, and 128. The L 2 
norm and H 1 semi-norm of the error on i.e., no MMPDE iteration, and F 5 \ 

i.e., 5 MMPDE iterations. The asymptotic order of the error is computed using two 
meshes with consecutive number of cells with respect to jt . 
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L 2 norm of the error 



H 1 semi-norm of the error 



Fig. 4.12: Example 2, mesh with N x N cells, for N = 8,16,32,64,128. Optimized for 
H 1 semi-norm instead of for L 2 norm. The L 2 norm and H 1 semi-norm of the error 
u — Uh after different numbers of MMP DE o uter iterations are plotted as functions of 
N. H 1 semi-norm based metric tensor (4.11 is used. 


Table 4.3: Example 2, mesh with N x N cells, for N = 8, 16, 32, 64, and 128. 
Optimized for H 1 semi-norm instead of for L 2 norm. The L 2 norm and H 1 semi¬ 
norm of the error on T^°\ i.e., no MMPDE iteration, and T^ 5 \ i.e., 5 MMPDE 
iterations. The asymptotic order of the error is computed using two meshes with 
consecutive number of cells with respect to H 1 semi-norm based metric tensor 

(4.1) is used. 



On mesh T ^ 
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